StatQuest生物统计学专题 - library normalization进阶之DESeq2的标准化方法
RPKM,FPKM,TPM等标准化方法还有那些问题？DESeq2的标准化方法的原理就是提高中等表达基因的地位一个例子
在上一节的StatQuest生物统计学专题中，我们简单直白的讨论了RPKM,FPKM,TPM的定义和生物学意义，明白了RPKM,FPKM,TPM标准化方法就是为了去除基因长度和测序深度对测序Read数的影响，见StatQuest生物统计学专题 - RPKM,FPKM,TPM，并且在一般意义上，更推荐使用TPM标准化方法。
然而不得不说明的是，诸如DESeq2、edgeR等差异表达分析软件都不是使用的RPKM,FPKM,TPM方法，为什么呢？
RPKM,FPKM,TPM等标准化方法还有那些问题？
其实TPM之类的标准化方法虽然解决了基因长度和测序深度的影响的问题，但还是不能解决一个问题：那就是测序文库组成不同造成的差异。
什么意思呢？
我们知道RPKM,FPKM,TPM是可以解决由于测序深度的差异而引起的Read数变异，如下图所示，样本2#的总Read数是样本1#的2倍，每个基因的Read数也是1#的2倍。我们知道这种Read数差异并不是基因表达的不同，而是由于测序深度不同所致，只需要将样本1#、2#除以各自的总Read数，那么这个Read数变异就会得到修正。
让我们再考虑两种新的情况，我们知道RNA-seq或其他高通量方法，往往是不同组织间的基因表达比较，比如肝细胞或脾细胞，他们会有一些各自特异性表达的基因，或者对于一些敲减基因的实验，有些基因在其中一个样本中高表达，而在另一个样本中不表达。
如下图所示的两个样本，两者的测序深度是一样的，总Read数相同。但是由于A2M基因在样本2#中不表达，导致563个Read会分配到样本2#的其他5个基因中去：如样本2#的基因A1BG的Read数是235，大于样本1#中的30。然而实际上这种差异并不是生物学效应所致，而是由于基因A2M被敲减所致，并不是A1BG等5个基因在样本2#中表达增加了，而是基因A2M在样本2#中表达减少了。
在这种情况下，这种测序文库组成不同的差异是RPKM,FPKM,TPM等方法无法解决的。
DESeq2的标准化方法的原理就是提高中等表达基因的地位
那么如何解决这个问题呢？
本周先看一下DESeq2是如何进行library normalization的，DESeq2的标准化方法共有7步，看起来很繁琐，但是原理很简单，它有一个贯穿始终的基本思想——提高中等表达基因的地位。
而且这7步只是为了得到一个标准化因子，并进行变换。
首先以下述数据集为例，共有3个样本，每个样本有3个基因：
第一步 对Read矩阵取对数变换
DESeq2默认是使用自然对数，也可以使用log2或log10。
第二歩 取各基因的平均数
要说明的是，由于样本1的基因1的Read数是0，所以在取对数时，它的值是-Inf（负无穷大），因此对基因1取平均数时就直接得出是-Inf即可。
为何要取对数？
其实是为了减少高表达异常值对标准化的影响，需要注意的是异常值不代表是错误值，只是说它相比数据趋势比较异常。
以原始表达矩阵为例，基因3的3个Read数分别是33、55、200，尤其是200，相比较整个表达矩阵来说是一个高表达异常值。
如果对原始表达矩阵求基因均值，那么基因3的均值是
(33+55+200)/3 = 96，而对于对数变换后的表达矩阵来说，基因3的均值是(3.5+4.0+5.3)/3 = 4.3； e^4.3 =73.7，73.7<96，也就是说对数变换后基因3的均数受到200的影响更小（这种取对数求得的平均数是几何平均数）。
第三步 过滤掉-Inf基因
将存在-Inf值的基因过滤掉，过滤掉的基因不再参与标准化因子的计算。
实际上，这一步是把在一个或多个样本中存在零表达的基因剔除。假定本实验是在比较不同组织细胞如肝细胞和脾细胞的表达量差异，那么这一步会剔除掉组织特异性表达的基因，而只保留管家基因——在不同细胞中都或多或少会表达的基因。
第四步 将对数矩阵减去对数均值，得到对数比值矩阵
将对数矩阵的每个Read分别减去此基因的对数矩阵。
对数相减的意义何在？
对数相减其实是真数相除，也就是：
log(reads for gene X)-log(average for gene X) = log(reads for gene X/average for gene X)注意：此时的average for gene X是几何平均数。
第五步 计算每个样本的对数比值矩阵的中位数
取中位数，而不是均值，也是为了进一步降低异常值的影响，具有较大表达差异的基因对中位数的影响甚微。考虑到绝大部分情况下，表达差异大的基因都是很少的，所以这个“中位数”更能代表的是中等表达基因或管家基因的情况。
第六步 将对数中位数转换为其相应的真数，得到各个样本的标准化因子
第七步 将原始表达矩阵除以这个标准化因子
原始矩阵的每个样本的全部Read数均除以各自的标准化因子（包括较大表达差异的基因）。
我们可以看到标准化之后，对于基因3来说，样本1#的Read数提升，而样本3#的Read数下降，基因3在3个样本中的表达其实是接近的。
做一下总结：
对数变换可以减少只表达在某些样本中的基因的影响，同时还可以减少异常值的影响（几何平均数）；
将在某些样本中表达量为0的基因剔除，不参与中位数计算，可以剔除特异性表达基因的影响；
中位数算法可以进一步降低高差异表达基因的影响，而提高中等表达的基因的地位；
标准化因子的生物学意义
其实这个标准化因子算法就是选出一个有代表性的gene X（其实是每个样本一个代表性gene X），而这个gene X的reads for gene X/average for gene X比值就是标准化因子。
只不过选取gene X的时候，通过对数变换和中位数的方法，更多的参考了中等表达基因和管家基因的数据趋势，而剔除了特异性表达基因和高差异表达基因的影响。
相比较RPKM,FPKM,TPM标准化方法是除以总Read数，DESeq2标准化方法是除以一个有代表性基因的Read数，只不过这个Read数进行了变换（它除以了几何平均Read数， reads for gene X/average for gene X）。因为更能处理存在特异性表达基因和高差异表达基因的数据。
一个例子
我们按照这个DESeq2的标准化方法的思想，对图2中的数据进行一个简单的标准化，没有完全按照上述的7步法，只是体会这种标准化的意思（结果没有大的差异，但是算法并不完全正确）。
#Gene Sample#1 Saple#2
A1BG 30 235
A1BG-AS1 24 188
A1CF 0 0
A2M 563 0
A2M-AS1 5 39
A2ML1 13 102
首先找到参照基因Gene X
对于Sample 1#来说，A1BG-AS1和A2ML1是中位数基因（不计算A1CF、A2M），而样本 2#也是同样。
为了简便，就用A2ML1基因了，其实结果是类似的。
原本应该根据ln(reads for gene X/average for gene X)的值计算中位数，这里直接根据原始Read值计算，会有一定的差异。
求解reads for gene X/average for gene X比值
由于，
average for gene A2ML1= (ln13+ln102)/2 = 3.12于是，
reads for gene A2ML1/average for gene A2ML1= (13,102)/3.12 = (4.16,32.66)也就是说两个样本的标准化因子分别是4.16和32.66。
进行标准化变换
将原始矩阵按照样本的不同除以各自的标准化因子，得下表：
可以发现，除了A2M基因有表达差异外，其他基因表达无明显差异。
#Gene Sample#1 Saple#2
A1BG 7.205869671 7.194095374
A1BG-AS1 5.764695737 5.755276299
A1CF 0 0
A2M 135.2301542 0
A2M-AS1 1.200978278 1.1939137
A2ML1 3.122543524 3.122543524
参考资料
StatQuest课程：https://statquest.org/video-index/
猜你喜欢
生信菜鸟团-专题学习目录（6）
生信菜鸟团-专题学习目录（7）
还有更多文章，请移步公众号阅读
▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。
▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。